function out = plognpdf(X, alpha, mu, sigma)
%   See "The Double Pareto-Lognormal Distribution - 
%   A New Parametric Model for Size Distributions", W.J Reed & M. Jorgensen,
%   http://citeseer.ist.psu.edu/691544.html
%
%   Created by: J. Huntley, 02/20/08
% 
%   Edited by Ben Lockwood, 03/26/2013
%   Changed to single Pareto log-normal. Limit of equation (23) from the paper as beta
%   goes to infinity.
%   Edited by Uwe Thuemmel, 2020, correction at extreme values. 
    
    X(X<0) = 0;     % replace negative elements with zeros
    arg2 = (log(X) - mu - alpha.*sigma.^2)./sigma;
    A = exp(alpha.*mu + 0.5*alpha.^2.*sigma.^2);
    out = alpha.*X.^(-alpha-1).*A.*normcdf(arg2); % equation (20)

    % correction at extreme values
    out(X==0) = 0;
    out(isnan(out)) = 0;
    out(out==Inf) = 0;
    out(out==-Inf) = 0;
end